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The Extended Parabolic Equation Method and Implication of 
Results for Atmospheric Millimeter- Wave 
and Optical Propagation 


1. Introduction 

During the three decades since the introduction of the stochastic parabolic wave equation by 
Tatarskii and Klyatskin [ Tatarskii , 1969], [Klyatskin, 1970], this theory of wave propagation 
through random media has found much success in the applications, from laser beam propagation 
through the atmosphere to image propagation. The applicability of this theory, however, is 
limited by the fact that it only describes situations in which the paraxial approximation prevails. 
Thus, the parabolic wave equation can be applied only when the smallest size of the permittivity 
fluctuations l 0 and the wavelength A satisfy the condition l 0 » A, i.e., the wavelength must be 
the smallest spatial scale in the problem. Although extensions of this theory have been advanced 
(see, for example, [ Saichev , 1980a], [Ostashev and Tatarskii, 1995]), none have been 
analytically treated at the level where corrections to the paraxial approximation have been 
identified and compared to paraxial results. 

The purpose of this memorandum is to provide a theoretical foundation for the extension of 
the theory which results in an ‘extended’ parabolic equation capable of treating wide-angle 
propagation in situations where the condition l 0 » A is violated, thus potentially making the 
extended theory formally applicable to millimeter wavelengths. Such a formulation is 
potentially required in propagation situations in which the operating wavelength is on the order 
of the spatial size of the inhomogeneities which make up the random medium. Such a theory 
will find application in the study of, e.g., millimeter wave propagation through atmospheric 
turbulence, in which the size spectrum of the turbulent eddys subtends a few millimeters, or in 
optical propagation through dust and sand storms where the wavelength can be on the order of 
the size of the inclusion. Most importantly, however, is the fact that solutions to the extended 
theory will provide a measure of correctness of the paraxial theory. 

In section 2, a first-principles derivation of the complete equations for forward and back 
scatter propagation are obtained from the stochastic Helmholtz equation, the various forms of 
which yield known results used in other such studies. Section 3 presents the derivation of an 
’extended’ stochastic parabolic equation from the general equations for forward scatter presented 
in Section 2. Here, the extended parabolic equation obtained is an operator equation which 
transcends the restrictions of the paraxial approximation but reduces to the latter in the 
approximation / 0 » A . 

Since the electric field as described by the extended parabolic equation is a random quantity, 
one can only deal with the related statistical quantities such as the first and second order 
moments. This is done in Section 4 where a general operator method is developed to treat 
generalized statistical moments within the extended parabolic equation. Solutions of the specific 
equations for the first order moment and second order moment (i.e., the mutual coherence 
function) is dealt with is Section 5 for the case of the random atmospheric permittivity field as 
described by the Kolmogorov spectrum of fluctuations. An appendix is provided which makes 
contact with the developments of Section 2 and other operator representations employed in the 
treatment of wave propagation problems. 


NAS A/TM— 2004-212944 


1 



2. Derivation of the Equations for Forward and Back Scatter 
Propagation form the Stochastic Scalar Helmholtz Equation 

Consider the propagation of a scalar electric field component of an electromagnetic wave in a 
random medium characterized by a stochastic permittivity e(r) = l-t-e(f) where e(r) is the 
random part. The Helmholtz equation governing the resulting stochastic electric field in the 
scalar approximation (i.e., neglecting wave depolarization) is then given by 

d + V 2 p E i x ’P) + k 2 E(x,p) = k 2 e(x,p)E(x,p ) (1) 

where the direction of wave propagation along the otherwise arbitrary x-axis is separated out 
from the three dimensional coordinate r = (x,p). The total stochastic field E(x,p ) can be 

decomposed into a forward propagating wave field E + (x,p ) and a backward propagating wave 
field E~(x,p), i.e., 


E(x,p) = E + (x,p) + E (. x,p ) (2) 

Similarly, one has for the derivative of the fields 

dE(x,p) _ dE + (x,p) | dE~(x,p ) 
dx dx dx 

As shown in Appendix 1, one can write these wave field components as an expansion into 
inhomogeneous plane waves, 


E± {x,p)=j JV (q)exp\iq ■ p±i(k 2 -q 2 y x\d 2 q 


( 4 ) 


Thus, 


~p~ = J 1 ^{^[^{k 2 -q 2 f ^] ex P iq ■ p ±i(k 2 - q 2 f~ x\d 2 q (5) 


/ 2 2\ 1/2 

This expression can be rewritten by expanding the factor (k —q~) 
series, 


in the integrand into a 
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— = J \ e + (cf) ±i\k-^q 1 -\ — 1 exp iq ■ p± i[k 2 —q xj clq 

= J J e± (q) ±{* 2 -^Vp + -" jexp iq- p± i(k 2 - ^ 2 ) 1/2 xjrf 2 ^ 

= J J e± (^)[± ? (^ 2 — Vp)' " exp iq- p±i(k 2 -q 2 f~x\d 2 q 

= ±i(k 2 — J j e ± (q)exp iq ■ p±i(k 2 -q 2 y xj d 2 q 

= ±i(k 2 -V 2 P ) 1/2 E ± (x,p) 

Therefore, applying this differential operator expression to eq. (3) gives 

d ~^~ = z t ^ 2 + V p)' [ E+ ( X 'P) ~ E ( X ’P)\ 

Differentiating this expression one more time, 

d 2 E(x,p) _ ./, 2 V 2 \'/ 2 dE + (x,p) dE(x,p) 
dx 2 l[l p) [ ax ax 

Making the identification [5 = [k 2 + ' and equating eqs. (3) and (7) yields 

^^ + ^^ = i P[E , {xrp) _ E -( X ,p)} 

Solving eq. (9) for dE + {x,p)/ax and substituting the result into eq. (8) gives 

Similarly, solving eq. (9) for dE~(x,p)/dx and substituting the result into eq. (8), 


d 2 E(x,p) 

dx 2 


f3 2 [E + (x,p)- E (x,p)\ + 2i(3 ^ E ^ 


Finally, substituting eq. (10) into eq. (2) and using eq. (1) gives 


2 i — -2 [3E (x,p) = -k 2 f5 l e(x,p)[E + (x,p)+ E ( x,p )] 


( 6 ) 


(7) 


( 8 ) 


(9) 


( 10 ) 


( 11 ) 


( 12 ) 
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Similarly, substituting eq. (1 1) into eq. (2) and using eq. (1), 

2 i -^-—^ + 2 [3E + (x,p) = k 2 /3~ l £(x,p)[E + (x,p) + E~(x,p)\ (13) 

C/yi 

Thus, eqs. (12) and (13) are, upon reinstating the expression for /), 

2 i — + 2[k 2 + Vp) V “£' + (x,p) = k 2 (k 2 + V 2 ) V ~e(x,p)[E + (x,p ) + E(x,p)\ (14) 

and 


2 i — -2 (k 2 + V 2 p f~E (x,p) = -k 2 {k 2 + V 2 p ) y "e(x,p)[E + (x,p) + E (x,p)\ (15) 

in which [k 2 + V 2 ^ is a differential operator and (/r + V 2 ) is an integral operator. (See 

Appendix 1 for these identifications.) If eqs. (14) and (15) are added together, one obtains the 
additional relation 


2* -~dy ^ +2 {k 2 + V 2 p f~\E + (x,p)-E (x,p)] = 0 (16) 

Equations (14) and (15) can also be written in terms of the Fourier transform of the field. 
Defining 


E(x,k ) 


1 


J £(x,p)exp(-/ic- p)d 2 p 


(17) 


Eqs. (14) and (15) become 


2 i ^ + 2 (k 2 - k 2 )' /2 E + (x, k ) = k 2 (k 2 - k~ 2 ) 1/2 e(x,ic)* [e + (x,k)+ E (x,/f)] (18) 


and 


2 i — -2 {k 2 -K 2 ^ 12 E (x,K) = -k 2 (k 2 -K 2 y ,2 £(x,K)*[E + (x,K)+ E (x,ic)] (19) 


where 

e(x,k)*E (x,k)=J \e(x,K r )E (x,k--7c')c/V (20) 
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denotes the convolution. Equations ( 14)-( 16), or (18) and (19) provide the foundations of many 
wave propagation studies [ Frankenthal and Beran, 1997], [Saichev, 1980a, b], [Malakhov and 
Saichev, 1980] and [Ostashev and Tatarskii, 1978,1979]. 


3. An Extended Stochastic Parabolic Equation 
for Wide-Angle Forward Propagation 

Forward scatter propagation in which the backscattered component is negligible, i.e., 
E~(x,p) « E + (x,p) is described by eq.( 14) in which E~{x,p) = 0, viz., 

2 i d E + l(k 2 + V 2 p f‘'E + (x,p) = k 2 (k 2 + V 2 ) V ~e(x,p)E + (x,p ) (21) 


This relationship, although capable of describing propagation situations in which the wave is 
scattered at angles up to n [ 2 with respect to the propagation direction, is still, strictly speaking, 
limited in its application to cases in which X < l 0 as noted in an analysis by [Ostashev and 
Tatarskii, 1995]. Thus, in this extension of the theory, the restriction X« l 0 is lifted only to be 
replaced by the less stringent condition X < / 0 . However, in what is to follow, situation in which 
A>/ 0 will also be considered as approximations. Equation (21) can be put into the form 
isomorphic to the standard paraxial parabolic equation by multiplying by 2 ik , rearranging some 
factors, and identifying a new set of differential and integral operators given by 



( 22 ) 


which allow eq. (21) to be written as 

dE + (x,p ) ~ r 

2 ik + U p E + (x,p ) - 2 k 4 V p {e{x,p)E + (x,p)} = 0 (23) 

Equation (23) can be called a generalized or extended stochastic parabolic equation in the 
differential operator U and the integral operator V . Reduction of this general relation to the 
paraxial parabolic equation is accomplished by expanding the operators U ~2k +V p and 

V p ~ l/2k 2 . Equation (23) then becomes 
dE + (x p) 

2 ik E7E1 + 2 k 2 E + (x,p) + W 2 E + (x,p) - k 2 e(x,p)E + (x,p) = 0 (24) 

ox 

Letting E + (x,p) = W (x,p)cxp(ikx) and substituting into eq. (24) gives 
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2 ik i—^l + V 2 W (x,p)-k%x,p)W ( x,p ) = 0 

ox 


( 25 ) 


which is the classical stochastic parabolic equation in the paraxial approximation. 

The extended stochastic parabolic equation eq. (23), just as its specialized counterpart, 
eq. (25), can be written as an equation governing the statistical moments of the electric field; this 
will form the subject of the next section. 


4. The Generalized Statistical Moments of the Electric Field 
Within the Extended Stochastic Parabolic Equation 

A. Derivation of the Fundamental Equations 

Consider the n + m th statistical moment of the electric field defined by 


/ n n+m \ 

Kn{^p v p 2 ,pp,p n+1 ---,p n+ J = ljjE(x,p i )YlE\x,Pj)) (26) 

\ <'=1 j=n + 1 / 

in which each of the transverse coordinates p i exists in a plane located at a distance x from the 
origin, at which a source is placed, in an otherwise arbitrary coordinate system. Defining the 
random operator 


+C. + 2 k\e(x„p,) 


Eq.(23) applied to E(x,p) becomes 


A,3,£(*>A) = ° 


Similarly, the complex conjugate of this relationship taken at (x,p\ gives 

d:, p e'(x.p^o 

Multiplying each of these equations by the remaining n + m - 1 th E fields gives 

n n+m 

n E (+pi){ D ,,p t E ( x ’pk)} n e*(x,pj)=o 


i=i 

tek 


j=n + 1 


(27) 


(28) 


(29) 


(30) 


and 
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( 31 ) 


n n+m 

IT E(x,pi) II E*(x,p j ){D* xp E*(x,p l j} = 0 

1=1 j=n + 1 

Summing eq. (30) over the n possible values of k, and summing eq. (31) over the m possible 
values of / and subtracting the second result from the first gives 


k=l 


n E (x,p){D xrp E(x,p k )}Y[ E *( x ’Pj) 


i= 1 
i*k 


j=n + 1 


l=l+n 


n n+m 

n E{x,p) n E *( x ’pj){ D U E *( x ’pi )} 


*=i 


j=n+l 


■ 0 (32) 


Using eq. (27) and rearranging operators as they operate only on E -fields with specific (x,p) 
coordinates gives, after some algebraic manipulations, 


2 ik — 

dx 


n n+m 

n em) n E\x,pj) 


;=i 


j=n + 1 


+ 


n n+m 

+ Up,, + ^V,e(x,p t ))- £((/;, + 2 k>v;e-(x,p,)) 

-j = 1 * 

n n+m 

n ^(+P,)n £ ‘(+p;) =o 


/=l+n 


i= 1 


7=«+l 


(33) 


Finally, defining the operator 


U s2 *£+£(U + 2 *+ + 2 k‘v; i e(x,p l )) 


1=1+ n 


Eq.(33) becomes 


where 


^nm & nm ’ Pn’ Pn + 1’ ’ Pn+m ) ^ 


n n -riri 

Snm ( + Pl’P2’' * ' ’Ph ’P 11 + 1>‘ ' ' 'Pn+m) = 1 1 Pf^'P/) J. J_ E (-+P/) 


1=1 


7=ii+l 


(34) 


(35) 


(36) 


NAS A/TM— 2004-212944 


7 



Equation (35) is an operator equation governing the product g nm of the n + m field values within 
the random medium. It is desired to obtain the ensemble average of the product 
(Snm (' ' ')) = T'nm (' ‘ ') given by eq. (26). Equation (35) cannot simply be ensemble averaged since, 

as the operator L nm is itself a random function of e{x,p) of which the values of the E -field are 
coupled, one has that (L nm g nn }j + (L m „^T nm . In order to obtain a closed equation for the field 

moment T nm , a stochastic operator method can be applied to eq. (35) which will yield an 
equation for T nm . 


B. Development of an Operator Method for the Equation Governing the Field Moments 

It is expedient to adopt the methods of [ Tatarskii and Gertsenshtein, 1963], and [Manning, 
1989] and decompose the operator L nm into its average and random parts, i.e., 



where 


and 



(38) 


L = 2k‘ 


n /?+/; ? 

'Lv P] e(x,p j )-J,Vp i e(x,p l ) 


. 7=1 


1=1+ n 



(39) 


Hence, eq. (35) becomes 


(( )<?,»» — 0 - ( 4 °) 
Ensemble averaging this relation yields 

(Am ) r „m + {KmSnm ) = 0 ( 4 1 ) 


Similarly writing 


8 , 


= r 

nm nm 




(8nm) = ° 


and substituting into eq. (41) gives 



r +(/. 


8 „ 


nm o nm 


= 0 . 


(42) 


(43) 
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Remembering that it is the goal of this development to obtain an expression for the general field 
moment T nm , one follows the development given in [Manning, 1989] and subtracts eq. (43) from 
eq. (35) and using eq. (42) once again obtains 


^nm^nm nm \~*nm nm \^nm&nm / ^ * 


(44) 


Combining the first and third members of this equation using the fact that 


L - L 


nm \ nm 


r = L r gives 

nm nm nm © 


$ n 


+ l r 

nm nm 



(45) 


One must now isolate the random quantity g nm by defining an operator L n ' m inverse to L nm , i.e., 
L~' m L nm = 1. Thus, operating on eq. (45) with L' vn yields 


Sn 


/V . ~ 

+ l:.l 


_ r 

nm nm nm 


-t 1 (L g ) = 0. 

nm \ nm<->nm / 


(46) 


Finally, operating on this relation with L nm , ensemble averaging, and solving the resulting 
expression for (L n J nm ) gives 



L L 1 

nm nm 


L L 1 L >r . 

nm nm nm / nm 


(47) 


Substituting this result back into eq. (43), one obtains for the equation governing T nm 


1 / \ 

1 -(L L 1 ) 

Vl L l L 

(L )- 

\ nm I 

\ nm nm j 

\ nm nm nm j 1 


(48) 


The solution of this operator equation, using eqs. (22), (38), and (39), gives an exact solution for 
the arbitrary field moments for wide-angle propagation through a random medium characterized 
by the stochastic pennittivity e{x,p). 

The general relation given by eq. (48) can be reduced to the parabolic equation for the field 
moments in the paraxial approximation in the case where X «/. In this case, one has the 
approximations for the operators U p ~ 2 k 2 + and V p ~ \/2k 2 used earlier. Equations (38) and 

(39) then become 





n+m 



l—l+n 


(49) 


and 
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( 50 ) 


L =k 

nm 


n n+m 

+J j e(x,p j )-J j e{x,p l ) 


L j = i 


1=1+ n 


One then employs two approximations; the first is given by 


L)+L 


nm / nm 


-1 


(51) 


where 



n n+m 

7=1 1=1+ n 

2 ik — 
dx_ 


— ]dx 

2 ik i 


and the second is given by 


1 


L L 1 

nm nm 


-1 


1 



(52) 


(53) 


under the assumption that the pennittivity fluctuations which enter eq. (50) are characterized by 
a zero mean value, i.e., (e^x,pj^ = 0. Equation (48) then becomes 


i2,4 + iv= 


Pi 


Zv 


1=1+ n 


L L L 

nm nm nm 


r =o, 

n nm. ‘ 


(54) 


where 


L L L ) « — 

nm nm nm / 


0 L /=! k= 1 


n n n n+m 

XX (e(x,p,)e(x',p())-X X {e(x,Pi)e {x',p'))- 

i= 1 /=l+n 

n+m n n+m n+m 

£ £ (e( x ,pj)p{x,p' k ))+ £ £ (e* [x,p j)e* (x',p') 


j=l+n k= 1 


7=l+n /=l+n 


dh', (55) 


which is the well known classical paraxial form for the problem [Tatarskii, 1971], [ Manning , 
1993]. It is interesting to note that the ‘geometrical optics’ approximation made in eq. (52) leads 
to the classical parabolic equation. Thus, one can envision a substantial extension of this 
development beyond that of the classical treatment if one is to use the entire form of the operator 

[L m )j , i.e., use the inverse of the operator (L nm j as solved in the paraxial approximation rather 
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than the geometrical optics approximation which was used above. This will form the subject of 
the next section. 

5. Solutions for the First- and Second-Order Statistical Moments: 

The Average Field and the Mutual Coherence Function 
for Wide-Angle Scattering 

A. First-Order Moment: The Average Field 

The first order moment of the random electric field in a plane transverse to the direction of 
propagation is defined through eq. (26) to be given by 

r H>(*; a) = ( E ( x 'A)) (56) 

which is a solution of the operator relation of eq. (48), viz., 

{(L, 0 )-[ 1 -(AA-;)((^o)}r 1 .=° (57) 

where, from eqs. (38) and (39), 

(L U) ) = 2ik^- + U p , 4 = 2 k 4 V p e{x,p), p = p r (58) 

The ability to proceed in an analytical fashion is dependent upon some simplifying 
approximations. Since L l0 = (L { 0 ) + L 10 , one has that 

<59) 

which allows one to approximately write 

1 — (AoAo) ~ 1 — (Ao)(Ao) (^0) 

where, analogous to the case met with in eq. (53), one has (V p e(x,p)^ ~ 0 . Using the results of 
eqs. ( 59) and (60), eq. (57) becomes 

{(i,o)-(i,o(Ao) T»)j r io =° ( 61 ) 

Employing the appropriate definitions of the operators, eq. (61) gives 
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2 ik ^ + U p -I (2 k\ e(x,p))\ 2 ik y+U p \ (2 k 4 V p e(x,p )) ) |r i0 (x,p) = 0 


This differential equation in the operators U p and V p must now be solved for the first-order 
moment r i0 (v,p). 

To this end, one must first deal with the factor 


2ik — + U p \ = G(x,p) 


which is the Green function of the operators 2ikd/dx + U , operationally defined by 


V P 


2ik — +2k 1+— f I rG(x,p) = <5(.r- x')d(p-p') 


? 1—7 ? 

where the definition of U p is used. Applying the approximation U ~ 2 k + V~ as well as the 
Fourier Transfonn relationship 


^ j ] G(x,p)exp(-ik p)d 2 p, 


eq. (64) can be written as 


dg \ 1 (2k 2 k 2 )c-( 1 p(*-v / )exp(-/A p) 

dx 2 ik ' ’ \2n) 2 ik 


The solution to this equation is 




exp -i(2.k 2 -K^(x- x') /2k exp {-ikp') 
(2n)\2ik) 


Finally, applying the inverse transform 


G(x,p) = J g(*, /c)cxp(i/c- p)d 2 


yields, upon evaluating the resulting integral in plane polar coordinates, 
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G(x,p ) = — — — U J exp\-i(2k 2 - k 2 )(x - x') /2k\exp{iKpcosO)KdKd6 

\2n)\2ik) Q 0 L v ’ 11 

= — | — J exp\-i(2k 2 - k 2 )(x - x') /2k\] Jkp)kcIk 
\2K\2ik) 0 L v ’ 1 i 


( 


exp [-ik(x-x')\ 


exp -ik(p- p'y 2(x - x') 


x- x 


= G(x,p-,x',p') (68) 

Thus, the third term within the brackets of eq. (62) can be written 

(•■•}^( 2 t 4 k } ( 2 * 4 h%p)) / 

= (2* 4 ) 2 | J G(x,p,x,p')(V p i(x,p)V p i(x,f>'))d 1 p'dx' (69) 

0 -OO 

Proceeding further, one now must deal with the operator products 


V p e{x,p ) = 


( 1 ] 

{ v 2 Y 
i+^ 

{2k 2 } 

k 2 

l ) 


e(x,p) 


(70) 


Since e(x,p) is a random function, it can be represented in the form of a Fourier- Stieltjes integral 
[Manning, 1993], i.e., 

e(x,p) = J exp (jit - p)dZ(x,K) (71) 


in which the spectral amplitude dZ(x,p ) is endowed with the same statistical properties as is the 
random function e(x,p) as will be shown in what is to follow. Applying eq. (71) to eq. (70) 
results in the following development: 


V p e{x,p) 



exp (iK- p)dZ{x,K) 

• ] j exp(zk- p)dZ[x,K ) 


- + ■ 


exp(zk- pyiZ(x,K) 


exp(zk- p]dZ(x,K) 


(72) 


NAS A/TM— 2004-212944 


13 



Thus, the ensemble averaged product appearing in right side of eq. (69) becomes, 


(V p £(x,p) v p ,£(x,p')) 



1- 


K 


y /2 ( 



Vl/2 

/ 


■exp(zk- p+iK' ■ p')^5Z(x,k)5Z(jc' 


( 73 ) 


One now makes use of the fact that the atmospheric permittivity fluctuation field e(x,p) is taken 
to be statistically homogeneous, characterized by a power spectral density <f> £ (x,k) in the 
transverse plane, and ^-correlated in the longitudinal direction; these circumstances allow one to 
write [ Manning , 1993] 

(dZ(x,K)dZ(x',ic^) = 2k8(x - x')8(k+ k')0 £ (x,k)5WV (74) 


Using this in eq. (73) and perfonning the integrations where possible yields 


(v p e(x.p)v p .i(x'.p)) 




exp[zk- p d ]<b E (x,K)d 2 K 


(75) 


where p d =p-p' is the difference coordinate. 

Equation (69) can now finally be evaluated by substituting into it eqs. (68) and (75); 
converting the integration in the p d -plane into one in plane polar coordinates and perfonning the 
associated angular integration gives 



K 


! Y 1 


exp [~ik(x - 



• J J 0 {Kp d ){x-xyUxp[-ikpl/2(x-x^]p d dp d ® e (x,K)dx'd 2 K (76) 

0 


Continuing on and perfonning the p d and x integrations gives 


( y-Y 
i-Y 

— V k ) 


\ f i ^ ^7 4 T \ 1 ^ / — \ 7 2 

Hdyr J 1 1- 7T ® e (x>*)d-K 


where the 5-function relation 


J <5(x - x'yix' = — 


(77) 


(78) 
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is employed. Finally, noting that the statistics governing the random field e(x,p) are not only 
homogeneous but also isotropic [ Manning , 1993], one has that O. (x, k) = O, (x, k) thus allowing 
the integral above to also be evaluated in plane polar coordinates and, performing the angular 
integration, yields 


(■■■) = -in 2 k 3 


( k iY 

V k , 


<& £ (x,k)kcIk 


(79) 


Returning to eq. (62) and, substituting eq. (79) into eq. (62) gives 


2ik — + 2k 2 
dx 


( v 2 ^ 2 

i+— f 

k 2 

v K J 


+ iK~k 


2 u 3 


( k 2 Y 
v 1 ; 


$> e (x,K)KdK 


r 10 (v,p) = 0 


(80) 


The general solution to this equation is unknown. However, for the plane -wave case, one has 
that 


1+- f 

k 2 

V K J 


( V 2 T 

r m (x,p)= i+tt r 10 (.r)=r 10 (v) 

1 k ) 


(81) 


since the plane wave will not possess any transverse variations. In this special case, eq. (80) 
becomes 


2ik — + 2k 2 + in 2 k 3 
dx 


( k iY 

l-^r 
\ k J 


<& £ (x,k)kcIk 


r ioM = ° 


the solution of which is 


TC ^ 

r i0 (jc) = r io (0)exp| ikx ~ — k 2 x 


( k 2 Y 
i-^ 
v k ) 


0 e (x,K)KdK\ 


(82) 


(83) 


In the limit where (l-K 2 /k 2 ^ ~1, eq. (83) becomes the well-known small-angle scattering 

result from the parabolic equation method [ Tatarskii , 1971], [ Manning , 1993]. 

The spatial spectrum of the atmospheric permittivity fluctuations is given by the Kolmogorov 
spectrum 


0 £ (k) = O.O33C £ V 11/3 , — < k < — (84) 

Y l 0 

where C 2 is the structure constant governing the strength or level of the turbulent fluctuations of 
the permittivity e, L 0 is the outer (i.e., largest) spatial scale of fluctuations and l 0 is the inner 
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(i.e., smallest) spatial scale. [ Tatarskii , 1971], [Manning, 1993] Typically in the open 
atmosphere, L 0 ~ 100 meters and / 0 ~ 1 millimeter. In order to avoid convergence problems with 
integrals such as that of eq. (83), and at the same time, incorporate the effects of the limits to the 
spectral interval, the Kolmogorov spectrum is augmented with a lower and upper ‘cutoff 
reflecting the finite outer and inner scales of turbulence. In this instance, one uses the associated 
modified von Karman spectrum 

4> e W = 0.033C e 2 (A: o 2 + r ! )'" ,6 exp(-^/^) (85) 


where the spatial frequency corresponding to the outer scale of turbulence is K 0 = In/L^ and that 
corresponding to the inner scale of turbulence is given by K m =5.92 // 0 [Manning, 1993], The 
use of this spectrum in eq. (83) results in an integration that is not analytically amenable. To this 
end, the following approximation procedure can be used [Ishimaru, 1977] 


0 £ (k-) = 0.033C £ 2 ( 

[K 2 + k 2 ) 11/6 exp (~K 2 /K 2 m ) 

= 0.033C; 

(K 2 + K 2 f U/6 - (7f 0 2 + k- 2 ) _I1/6 + (k 2 + K -2 )” 1 1/6 exp(-?c 2 /?c 2 ) 

= 0.033C; 

(K 2 + K 2 y U/6 - (K 2 + K 2 y nb {l-exp(-K 2 / K 2 m )}] 

-0.033C; 

\k 2 + K-y U/6 - K n ' 3 {l-OX V (-K 2 /K 2 m )}\ 


( 86 ) 


where the last line results from the fact that jl - exp (-K 2 lK 2 m )} has a non zero contribution only 
when k»K q . Using this representation in the integral indicated in eq. (83) and evaluating, 
using Mathematica [Wolfram, 1999], the resulting expression about the singularity at K = k and 
retaining the Cauchy principal value yields 


f 10 (x) = r io (0)cxpl ikx - C 2 x\ 


0.5439£ 1/3 (-l) 1/6 exp| 


( 7,2 


11 k 2 ^ 


V < 


K 


m J 


-0.08 14/c 1/3 77T exp 


( 7,2 3 


V K mJ 


+ 0.0911 k 2 K~ 5/3 


( 

1, 

l;I 



V 

6 

k2 ) 

IJ 


(87) 


where T(- ••,•••) is the incomplete gamma function. At the outset, this expression can be 
simplified by noting that one always has that k»K 0 . The hypergeometric function ,Ff ■ •) thus 
can be approximated by unity, the first term in its series expansion. Also, for the same reason, 
the second term within the braces is negligible with respect to the third term. One then has 


r, 0 (v) = r io (0)cxpl ikx - C 2 x 


0.5439£ 1/3 (-l) 1/6 exp| 



•2 \ 


K, 


+ 


m J 


+0.0977 k 2 K~ 5/2 }\ (88) 
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Finally, replacing the incomplete gamma function in terms of the slightly more manageable 
Kummer confluent hypergeometric function [ Grads hteyn and Ryzhik, 1980, eq. (8.351.4)], 



K 


Y/6 


k: 


exrt 




(89) 


yields 


r i0 (x) = r io (0)exp ikx - C 2 ;t{0.5439fc 4 kJ 1/3 U 


3 


17 _k*_ 

^ s ’ 2 

V 6 


+ 


+0.0977/c 2 J fif“ 5/3 }] (90) 


This expression will now be examined in the limits k > K m and k < K m . In the first case, the 
confluent hypergeometric function becomes [Abramowitz and Stegun, 1965, eq. (13.5.2)] 


ufi, 


17 

~ 6 ~ 


K, 


Y 1 


m J 


K~ 


> k > 


V 


In this case, eq. (90) becomes 

r i0 (x) = r io (0) exp[/fcc - 0.0911k 2 K~ 5/3 C 2 x + 0.5439& 2 k; 5/3 C £ 2 x], k > K m (91) 

This result agrees with that obtained from the small-angle scattering form of the theory; see, e.g., 
\Ishimaru, 1977] in which the second tenn of eq. (91) is neglected since L 0 » /„ for atmospheric 
turbulence. 

In the other limit in which k< K m , one now employs the approximation [. Abramowitz & 
Stegun, 1965, eq. (13.5.6)] 



r i7 

k 2 ^ 


"if 


u 

i — 


~ T 




1’6’ 

Kn; 


v 6 , 

if- 

K 7 


k < K.„ 


Substituting this result into eq. (90) yields, after simplification, 

r i0 (x) = r io (0)exp[ifcc-0.0977A: 2 C £ X 5/3 x- 0.443 lk 1/3 C 2 x - 

-0.2558 ik 1/3 C 2 x\, k<K m (92) 


The third and fourth terms within the exponential represent the wide-angle scattering corrections 
to the mean field. The imaginary fourth term describes the diffraction phenomena inherent in the 
scattering mechanism. These correction terms are only appreciable in the event that the 
wavelength X~L 0 . In the case of atmospheric turbulence, in which 7«L 0 , these correction 
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terms become negligible. This result is, unfortunately, not so much a testament to the extended 
theory as presented here as it is for the applicability of the small-angle scattering approach within 
the context of Kolmogorov turbulence spectra. In situations in which the outer scale of 
turbulence is on the order of the wavelength, these additional perturbations could easily be 
identified. 

In many instances within the literature, one finds use of the refractive index structure 
constant C] which is related to that of the permittivity C £ 2 by the relation C ; 2 = 4C 2 . In addition 
to this, these results also appear in the form of the coherent intensity I c (x) defined by 

/ c (.t) S |(£(x,0)>| ! =|r io (x)f (93) 

Thus, in terms of C 2 , the coherent intensity associated with eqs. (91) and (92) is given by 

I c (x) = |r io (0)| 2 exp[-0.78 16 k 2 C 2 E K^x + 4.35 12 k 2 C^; n 5/3 x], k > K m (94) 

and 

I c {x) = |r io (0)| 2 exp[-0.7816^ 2 CX‘ 5/3 v - 3.5448fc 1/3 Cpc - 

-2.0464 ik 1/3 C*x\, k<K m (95) 

B. Second-Order Moment: The Mutual Coherence Function 

The second order moment of the random electric field in a plane transverse to the direction of 
propagation is, from eq. (26), given by 

T n(- X ’PvP2) = ( E { x 'A) E X x '’pj) (96) 

which is a solution eq. (48), in this case given by 

|(A,)-[ 1 -(A,i,l)r'{A,4:fi)} r ,,= 0 

where from the definitions of eqs. (38) and (39), 

( 98 ) 

and 

L u =2k‘[\ (L„> = 0 (99) 


NAS A/TM— 2004-212944 


18 



As with the case for the first order moment, an assumption is made at the outset to render the 
problem analytically tractable. In particular, one has that 


i-{A4"!)] ' *[i-(A,}(A,)’ l J 1 =1 (ioo) 

which allows eq. (97) to become 

{(A,)-(A,f!A,)}r,,=o (ioi) 

which, upon using the operator definitions, is 

2-*|t+A -u;,-{(2^{v f i(x,p,)-r f t(x,p,)}) ■ 

^2ikj- + u fi -u;^ ( 2 r{V p £(x,A)-V i ;/(^ft)}))r'M =0 (102) 



Proceeding as in Section C.l and using the approximations U p ~2 k 2 + V 2 p for these operators, 
eq. (104) becomes 

| 2 ^Jc + V ^ ~ V 12\ G ( X ^P2) = S ( X - X 1 S {A-P' 1 ) S {P2-P , 2) ( 106 ) 

Defining the two-dimensional Fourier Transform g(x,K v K 2 ) of G(x,p ] ,p 1 ) by 
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g(x,K v K 2 ) = M- J ] G(x,p v p 2 )exp(-iK l p l -iK 2 -p 2 )d 1 p l d 2 p 2 
\2 k J ^ ^ 


the transform of eq. (106) is 

dg f 1 


dx l 2 ik 


2 2 \ ( 1 Y S(X-X') 

2 . k - exp(f/c, ■ p, + ik 2 ■ p 2 ) 


(108) 


and the solution of which is 


cxpi-if/c, 2 - kI)(x - x') ilk 

g{x,K v K 2 ) = . exp(i/c l ■ p; + iK 2 ■ p 2 ) 

[2k) (2 ik) 


(109) 


Applying the transform inverse to that of eq. (107) to eq. (109), converting to plane polar 
coordinates, and evaluating the associated angular integrals gives 

G'(^PpP 2 )=^J^JJ J exp[-/(x - x')Kf/2k]S 0 (K t p dl ) • 

•exp[-/(v - x')k 2 /2k ] )i Q [K 1 p cl2 )K x k 2 cIk { cIk 2 (110) 

where p di = |p ( . — p'|. Finally, performing the remaining integrals over the spatial frequencies 
yields for the Green function 

G(x,p v p 2 ) = - ex P W(Pi - P') - (Pi - Pi))/ 2 ( x - *0] 


G[x,p v p 2 ,x,p[,p' 2 ) (111) 


Therefore, the fourth term in eq. (102) is given by 


•) = ((2^{F p e(.r, a) - V;/ (x,p 2 )}){- ■ ■y 1 (2k 4 {v p e(x,p l ) - V p */ (x,p 2 )})) 

= ( 2kA ) \ J J G(x,p p p 2 ;x(p;,p;)(l/ p y p; e(x,A)e(x,p;) 


- v P y P 'A x ’p (*' &) - v P y P £ (^p 2 )e( x ' pd + 


+ V P 2 V l'f ( X ' P 2 ) g * (*’ P 2 ) )' d '~p[ d P'l dX ' 


( 112 ) 


As was done in Section (5A) above, one now uses the Fourier-Stieltjes representation for the 
random functions e(.r,p); thus, the result of eq. (72) is now employed in eq. (1 12). Additionally, 
in order to carry out the evaluation of the ensemble averages over the spectral amplitudes, one 
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assumes statistical homogeneity in the transverse plane and d-corrclatcdncss in the longitudinal 
direction. The result of implementing this procedure on eq. (112) gives 



■ [exp (i k-(p,- A')) - ex P ( z * ' ( A - Pi)) ~ 

-exp (iK- (p 2 -pj’)) + exp(iK’- (p 2 - p' 2 ))\d 2 Kd 2 p'd 2 p' 2 dx' (113) 


Substitution of eq. (Ill) into eq. (113) results in an integral which is unwieldy but nonetheless 
straightforward to evaluate and follows along the same lines as the treatment of eq. (76). Thus, 
integrating over the spatial variables and using eq. (78) gives 


(•••} = -ink 3 



k 2 Y 

T2 <& e (K)[l-CXp(-iK- (p 2 -p l ))\d 2 K 
k ; 


= -iln 2 e] 
0 


( 

1 

v 



° e W[ 1 - J oWP2-Pl|)]^*' 


(114) 


where statistical isotropy is assumed in arriving at the last result. 

Finally, substituting eq. (114) into eq. (102) gives for the differential equation governing the 
wide-angle mutual coherence function 


d 

f 

V P 1 

2ik — + 2k 2 

1 + 


dx 

V 

k 2 

/ 


-2 k : 


V 

l+-£ 


2 Y 2 


+ 


+i2n k 


2 7,3 


( 


1-^ 

v k Y 


V 1 


° e W[ 1 - J o( K |P 2 -Pi|)]^ K ' 


r(x;Pi,p 2 ) = 0 (115) 


The solution to this operator equation is not known. However, in the plane wave case in which 
there are no transverse variations, expressions analogous to eq. (81) prevail, i.e., 


f 

1 + 

V 





->1 


(116) 


thus allowing eq. (1 15) to become 


2ik— + i2n 2 k 3 
dx 


( 


i-4 

k 2 


V 1 


° e W[ 1 - J oWp 2 -PiH^*' 


nx'PvPi) ® 


(117) 


the solution of which is 
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(118) 


r(x;p d ) = r(0;p d )exri-n k x 


( 


1-^ 

k 2 


V 




The application of the Kolmogorov spectrum given by eq. (84) can be used in eq. (118) since 
the structure of the integrand is such that no convergence problems will arise due to the neglect 
of the spectral cutoffs inherent in the use of the modified von Karman spectrum, eq. (85). 
(Besides, the use of the complete modified von Karman spectrum in eq. (118) results in an 
integral that cannot be analytically evaluated.) Thus, substituting eq. (85) into eq. (118) and 
performing the required integration in Mathematica [ Wolfram , 1999] gives 


r ii(^A/) = r n( () -A/) cx n 


-(0.8861- 0.5 1 16 i)k V3 C 2 x( 1 - J 0 (kp d )) - 


( 

-0.0271 k'Chtfx ,F 2 1; 

V 


17 17. k 2 p ; f 
6 ’ 6 ’ 4 


(119) 


in which , F 2 {- ■ •) is a generalized hypergeometric function. Relating this generalized 
hypergeometric function to a more familiar function, one makes use of the relation [Gradshteyn 
and Ryzhik (1980), eq. (8.574.3)], 


ft 


17 17 k'-p'd 


1; 6 ' 6 ; 


7-1 



-11/3 


i k Pd) 


( 120 ) 


where s 8 , Q (- • •) is a Lommel function. Hence, eq. (119) becomes 


rn(^A,) = r n (0.p, / )exp[-(0.8861-0.5116i)*: ,,3 C;. l (l-J 0 (< : p,))- 

- 03M3k' l 'Clxs % fip d ) 


( 121 ) 


This expression will now be examined in the limits kp d > 1 and kp d <1. In the first limit, one 
employs an asymptotic result obtained from [ Gradshteyn and Ryzhik, 1980, eq. (8.570.2)], viz. 


s% 0 ( k P d ) “ ( k P d ) 


5/3 


F 

3 J 0 


( 5 5 4 

v ’ 6’ 6” k 2 p 2 

(kp d ) 5/3 


kp d > 1 


and, of course, the result 1-J Q (kp d ) =1 to obtain (neglecting terms on the order of k l/i since 
&<p d ) 


r n(*’P<i) = r n(0’Pjexp[-0.3643k 2 C;pf x ], k P d > 1 


( 122 ) 
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In tenns of the refractive index structure constant Cf t , this yields 


r„(*,A,) = r M (0, ft )exp[-1.4572*: 2 C, 2 pf x], kp d >1 (123) 

which is the familiar plane wave result from small-angle scattering theory [ Tatarskii , 1971], 
[Manning, 1993]. 

In the other extreme where kp d < 1, one has 


1 J o(kp d ) 


1 2„2 
k Pd 


as well as [Gradshteyn and Ryzhik, 1980, eq. (8.570.1)] 


s 8/ n ( k Pd)~ 77 (■ k Pd ) 


11/3 


11 


from which eq. (121) gives 

r„(*,A,)» r n (0.p„)exp[-0.2215i 7 ' J C]p]x + 0.1279 ik^p^x - 

-0.027 fit 4 C;p“ /3 x] (124) 

or, in tenns of C \ , 

r„(^P,,) = r i ,(0.p„)exp[-0.886* 7 ' 2 C i ;p> + 0.51 16i k^C]p]x - 

-0.108 4k 4 C;p l d y3 x] (125) 

Hence, as with the case of the first order moment of the wave field, the corrections to the second 
order moment due to the use of this extended theory are negligible in the case of atmospheric 
turbulence at optical frequencies. Only in the case of millimeter waves could these small 
conections be noted. 


6. Summary and Conclusions 

Beginning with the stochastic Helmholtz equation, relations are obtained which transcend 
the paraxial approximation as far as diffraction phenomena are concerned but are limited to 
situations prescribed by / 0 > X as far as propagation phenomena is concerned. With this situation 
notwithstanding, a resulting extended stochastic parabolic equation for wide angle propagation 
has been derived which goes beyond the paraxial approximation of the classical small angle 
equation and reduces to it in the appropriate limits; this equation is an operator equation which is 
isomorphic to the classical small angle scattering equation. After the development of an operator 
method to treat such an equation for the statistical moments of the wave field, solutions of first 
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and second moments were obtained in the case of atmospheric turbulence. It is found that the 
corrections made to the classical results are insignificant in the atmospheric case. 

The Kolmogorov spectrum which is used to represent atmospheric turbulence is such that 
the contribution of the spectral frequencies which approximately correspond to the wavelength is 
relatively small as compared to those at the smaller spatial frequencies (which correspond the 
outer scale of turbulence). This gives rise to the insignificant levels of the correction terms at 
nominal operating wavelengths. These results support the use of the paraxial approximation of 
the stochastic parabolic wave equation at optical and millimeter wavelengths propagating 
through atmospheric turbulence. This may not be the case for propagation through atmospheric 
aerosols which possess an entirely different spectral behavior [Manning, 1993]. This case should 
be considered in further works on this subject. 
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APPENDIX 1 

Plane Wave Expansions and Various Operator Representations 


The expansion of an outgoing spherical wave 


, exp (ikr) „ - 

go{ x ’P) = . r = xx + p 

4 nr 

which is a solution of the Helmholtz equation 

(V 2 + ^ 2 )^o(-TP)=0 

is given by 

8o( x >P) = J A{q)txp(iq ■ r)d 3 q 


(Al) 


(A2) 


(A3) 


the spectrum A(cj ) of which can be found by inverting eq. (A3), i.e., 

J g 0 (x,p)exp(-iq' ■ r)d 3 r = J J A(q)exp(iq ■ r-iq ■ r)d 3 qd 3 r 

= (2tt) 3 J A(q)8(q-q')d 3 q 
= (2nfA(q) 


(A4) 


The left side of eq. (A4) is easily found by converting it to an integral in spherical coordinates, 
viz., 


J go(*>P)exp(-/<f- r)d 3 r = jL ex P(^0 exp (-/<f • 


' 3 


Atit J 


)d 


r)a r 


=~i J I 


exp (ikr) 


An 


cxp(-/V/Vcos 0) r 2 sin 0 dr dOd(j) 


0 0 0 


= — } J exp (ikr - iq'ra)r dr da 

2-i o 

= — — J 2i exp(ikr) sm(q'r) r dr 
2iq o 

1 


/2 / 2 
q —k 


(A5) 
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Thus, from eq. (A4), 


A(q) = A(q ) 


-T 

r i ) 

v 2ttJ 

2 7 2 

[q -k J 


Thus, eq. (A3) can be written 


8o(X'P) 



(A6) 


(A7) 


From this relationship, many operator representations can be obtained which are used in this 
work as well in other works. 

To this end, writing eq. (A7) in Cartesian coordinates gives 


8o{x,p) 



ex P [i(q,x + q p -p)\ 


2 . 2 

C l, + C ln 


dq x d 2 q 


performing the q x -integration by closing the resulting contour integral in the complex plane in 
the x > 0 region yields 


8 0 {x,p) 


87T 


II 


9p<k 2 


exp 

_*■ — * / 2 2 \ 

iq p - p + ix(k -q-J 


k 2 ~q 2 p 



Thus, one has in tenns of the two-dimensional plane-wave spectrum 


8o(x,P)= JJ ko(<7 P )e x P 


iq t 


■ p+ ix[k 2 



(A8) 


(A9) 


where 



(A10) 


In the ‘transverse’ case where the spectrum is considered across the interval p-p 0 in the 
transverse plane at .r = 0 , eq. (A9) becomes 
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^o(O^p-Po) = JJ g 0 (q P )^p[iq p -(p-po)\d 2 q p 

=(^)ll ( k 2 -y 2 P ) cx pf% ■ {p-m\ 

= (^)ll [ A ' 2 +\%~- • • j ex p[% • (p - Po)\ d \ 

= (^t)U^ 2 + ^t) +••• ex Pfe ' (P~Po)\d 2 q p 

= ^^rjlJ ( k2 + V p) V "exp [iq p ■ (p-p 0 )\d 2 q p 

= (ii /) A2 + ' H exp [^ ‘ (P~pA dlq p 

=(y* a+v ;r^-^ o) 

Thus, one has the operator representation 

2ig 0 (0,p - ft,) = (<r + Vjfi^p-ft,) (All) 

Additionally, operating on eq. (All) with the transverse Helmholtz operator k 2 + V 2 gives 

2 i(k'- + Vj)« o (0 ,p-p 0 ) = (* 2 + Vlfs(p-p„) (A12) 

Finally, returning to eq. (A 10), one can consider the derivative 

= U 8o(q P ){i(k 2 -q 2 p ) V2 }eW iq p • p + ix(k 2 ~q 2 p ) V ~] d 2 q p 
= *'11 gofeX^-Vpf'exp iq p ■ p + ix(k 2 -^)'] d 2 q p 

= i (k 2 -V 2 p fg 0 (x,p-p 0 ) (A 13) 

Equations (A1 1)-(A13) fonn the basis of the theory developed by [Babkin and Klyatskin, 1980] 
and [Babkin, Klyatskin, and Lyubavin, 1980] to which the methods of invariant imbedding were 
applied. 

The operators (k 2 + and (P + V^) can be given analytical representations using the 
foregoing. In particular, consider an arbitrary function F(p) and the definition of the ^-function 
to fonn 
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(A14) 


(k 2 + vlf F{p)= J (k'- + vifs(p-py(p')d'-p' 
= J Kip-pytpyy 


where 


K(p-p') s (k'- + Vlfs(p-p') 

But using the result of eq. (A 12) in eq. (A 15) yields 

K(p-p') = 2i(k 2 +V;)g n (0.p-p') 

Thus, returning to eq. (A 14) and using the definition of eq. (Al), 

i(k > +K y m= ± J ( ^ + vy^-f) 




In _ 


P-P 


/ 2 o \ — 1/2 

which is a differential operator. Similarly, for the operator (& +V) , 

(k 2 +v 2 p y V2 F(p)= J (k 2 + v 2 p y 1,2 8(p-p')F(p')d 2 P ' 

= J K\p-p')F(p')d 2 p' 


(A 15) 


(A 16) 


(A 17) 


where 


K\p-P’) = (*r + Vlfs(p-p')=2ig 0 (0,p-p') (A 18) 

where the last result is by eq. (Al 1). Again, remembering the definition of eq. (Al), eq. (A17) 
becomes 


;(* J +vjr' ,2 F(p) 



txp(ik\p-p'\) 


P-P 


F{p')d 2 p' 


(A 19) 


which is an integral operator. 

Both eqs. (A 16) and (A 19) are used extensively by [Saichev, 1980a,b] and [Malakhov and 
Saichev, 1980]. 
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